!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine weno (ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz, x, y, z, dx, dy, dz, dt, nx, ny, nz,                              &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd,                                       &
   fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: dx, dy, dz, dt, x(nx), y(ny), z(nz)
   real(8), dimension(nx, ny, nz), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, po, gx, gy, gz
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, nz, n_neighs), intent(inout) ::                                                                          &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd
   real(8), dimension(ny, nz, n_neighs), intent(inout) ::                                                                          &
      fzro_bnd, fzmx_bnd, fzmy_bnd, fzmz_bnd, fzbx_bnd, fzby_bnd, fzbz_bnd, fzen_bnd, fzpo_bnd

   integer(4) :: im1, jm1, km1
   real(8), dimension(nx, ny, nz) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre
   real(8), dimension(nx, ny, nz) :: fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo
   real(8), dimension(nx, ny, nz) :: fyro, fymx, fymy, fymz, fybx, fyby, fybz, fyen, fypo
   real(8), dimension(nx, ny, nz) :: fzro, fzmx, fzmy, fzmz, fzbx, fzby, fzbz, fzen, fzpo
   real(8), dimension(nx, ny, nz) :: bx_resis, by_resis, bz_resis, en_resis

   real(8) :: ch, ch2, cp2
   real(8) :: dt_dx, dt_dy, dt_dz, dx2, dy2, dz2

!----------------------------------------------------------------------------------------------------------------------------------|
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   dz2 = dz * 2.0d0
   dt_dx = dt / dx
   dt_dy = dt / dy
   dt_dz = dt / dz
   ch = cfl * min (dx, dy, dz) / dt * div_clean_flag
   ch2 = ch * ch
   cp2 = -dt * ch2 / log (cd) + eps

!----------------------------------------------------------------------------------------------------------------------------------|
!  first step of 2nd order Runge Kutta method
!----------------------------------------------------------------------------------------------------------------------------------|
   call weno_lx (fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo,                                                             &
      ro, mx, my, mz, bx, by, bz, en, po, nx, ny, nz, ch2)
   call weno_ly (fyro, fymy, fymx, fymz, fyby, fybx, fybz, fyen, fypo,                                                             &
      ro, my, mx, mz, by, bx, bz, en, po, nx, ny, nz, ch2)
   call weno_lz (fzro, fzmz, fzmy, fzmx, fzbz, fzby, fzbx, fzen, fzpo,                                                             &
      ro, mz, my, mx, bz, by, bx, en, po, nx, ny, nz, ch2)

   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            ro_pre(i, j, k) = ro(i, j, k) - dt_dx * (fxro(i, j, k) - fxro(im1, j, k)) -                                            &
                                            dt_dy * (fyro(i, j, k) - fyro(i, jm1, k)) -                                            &
                                            dt_dz * (fzro(i, j, k) - fzro(i, j, km1))
            mx_pre(i, j, k) = mx(i, j, k) - dt_dx * (fxmx(i, j, k) - fxmx(im1, j, k)) -                                            &
                                            dt_dy * (fymx(i, j, k) - fymx(i, jm1, k)) -                                            &
                                            dt_dz * (fzmx(i, j, k) - fzmx(i, j, km1))
            my_pre(i, j, k) = my(i, j, k) - dt_dx * (fxmy(i, j, k) - fxmy(im1, j, k)) -                                            &
                                            dt_dy * (fymy(i, j, k) - fymy(i, jm1, k)) -                                            &
                                            dt_dz * (fzmy(i, j, k) - fzmy(i, j, km1))
            mz_pre(i, j, k) = mz(i, j, k) - dt_dx * (fxmz(i, j, k) - fxmz(im1, j, k)) -                                            &
                                            dt_dy * (fymz(i, j, k) - fymz(i, jm1, k)) -                                            &
                                            dt_dz * (fzmz(i, j, k) - fzmz(i, j, km1))
            bx_pre(i, j, k) = bx(i, j, k) - dt_dx * (fxbx(i, j, k) - fxbx(im1, j, k)) -                                            &
                                            dt_dy * (fybx(i, j, k) - fybx(i, jm1, k)) -                                            &
                                            dt_dz * (fzbx(i, j, k) - fzbx(i, j, km1))
            by_pre(i, j, k) = by(i, j, k) - dt_dx * (fxby(i, j, k) - fxby(im1, j, k)) -                                            &
                                            dt_dy * (fyby(i, j, k) - fyby(i, jm1, k)) -                                            &
                                            dt_dz * (fzby(i, j, k) - fzby(i, j, km1))
            bz_pre(i, j, k) = bz(i, j, k) - dt_dx * (fxbz(i, j, k) - fxbz(im1, j, k)) -                                            &
                                            dt_dy * (fybz(i, j, k) - fybz(i, jm1, k)) -                                            &
                                            dt_dz * (fzbz(i, j, k) - fzbz(i, j, km1))
            en_pre(i, j, k) = en(i, j, k) - dt_dx * (fxen(i, j, k) - fxen(im1, j, k)) -                                            &
                                            dt_dy * (fyen(i, j, k) - fyen(i, jm1, k)) -                                            &
                                            dt_dz * (fzen(i, j, k) - fzen(i, j, km1))
            po_pre(i, j, k) = po(i, j, k) - dt_dx * (fxpo(i, j, k) - fxpo(im1, j, k)) -                                            &
                                            dt_dy * (fypo(i, j, k) - fypo(i, jm1, k)) -                                            &
                                            dt_dz * (fzpo(i, j, k) - fzpo(i, j, km1))
         enddo
      enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis, ro, mx, my, mz, bx, by, bz, en, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term(mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre,                                                &
      ro, mx, my, mz, bx, by, bz, po, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,                                          &
      ch2, cp2, dx, dy, dz, dt, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value for fine-coarse block boundary
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
      fxro_bnd(j, k, 1) = fxro(ng, j, k) * 0.5d0
      fxro_bnd(j, k, 2) = fxro(nx - ng, j, k) * 0.5d0

      fxmx_bnd(j, k, 1) = fxmx(ng, j, k) * 0.5d0
      fxmx_bnd(j, k, 2) = fxmx(nx - ng, j, k) * 0.5d0

      fxmy_bnd(j, k, 1) = fxmy(ng, j, k) * 0.5d0
      fxmy_bnd(j, k, 2) = fxmy(nx - ng, j, k) * 0.5d0

      fxmz_bnd(j, k, 1) = fxmz(ng, j, k) * 0.5d0
      fxmz_bnd(j, k, 2) = fxmz(nx - ng, j, k) * 0.5d0

      fxbx_bnd(j, k, 1) = fxbx(ng, j, k) * 0.5d0
      fxbx_bnd(j, k, 2) = fxbx(nx - ng, j, k) * 0.5d0

      fxby_bnd(j, k, 1) = fxby(ng, j, k) * 0.5d0
      fxby_bnd(j, k, 2) = fxby(nx - ng, j, k) * 0.5d0

      fxbz_bnd(j, k, 1) = fxbz(ng, j, k) * 0.5d0
      fxbz_bnd(j, k, 2) = fxbz(nx - ng, j, k) * 0.5d0

      fxen_bnd(j, k, 1) = fxen(ng, j, k) * 0.5d0
      fxen_bnd(j, k, 2) = fxen(nx - ng, j, k) * 0.5d0

      fxpo_bnd(j, k, 1) = fxpo(ng, j, k) * 0.5d0
      fxpo_bnd(j, k, 2) = fxpo(nx - ng, j, k) * 0.5d0
   enddo
   enddo

   do k = 1, nz
   do i = 1, nx
      fyro_bnd(i, k, 1) = fyro(i, ng, k) * 0.5d0
      fyro_bnd(i, k, 2) = fyro(i, ny - ng, k) * 0.5d0

      fymx_bnd(i, k, 1) = fymx(i, ng, k) * 0.5d0
      fymx_bnd(i, k, 2) = fymx(i, ny - ng, k) * 0.5d0

      fymy_bnd(i, k, 1) = fymy(i, ng, k) * 0.5d0
      fymy_bnd(i, k, 2) = fymy(i, ny - ng, k) * 0.5d0

      fymz_bnd(i, k, 1) = fymz(i, ng, k) * 0.5d0
      fymz_bnd(i, k, 2) = fymz(i, ny - ng, k) * 0.5d0

      fybx_bnd(i, k, 1) = fybx(i, ng, k) * 0.5d0
      fybx_bnd(i, k, 2) = fybx(i, ny - ng, k) * 0.5d0

      fyby_bnd(i, k, 1) = fyby(i, ng, k) * 0.5d0
      fyby_bnd(i, k, 2) = fyby(i, ny - ng, k) * 0.5d0

      fybz_bnd(i, k, 1) = fybz(i, ng, k) * 0.5d0
      fybz_bnd(i, k, 2) = fybz(i, ny - ng, k) * 0.5d0

      fyen_bnd(i, k, 1) = fyen(i, ng, k) * 0.5d0
      fyen_bnd(i, k, 2) = fyen(i, ny - ng, k) * 0.5d0

      fypo_bnd(i, k, 1) = fypo(i, ng, k) * 0.5d0
      fypo_bnd(i, k, 2) = fypo(i, ny - ng, k) * 0.5d0
   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      fzro_bnd(i, j, 1) = fzro(i, j, ng) * 0.5d0
      fzro_bnd(i, j, 2) = fzro(i, j, nz - ng) * 0.5d0

      fzmx_bnd(i, j, 1) = fzmx(i, j, ng) * 0.5d0
      fzmx_bnd(i, j, 2) = fzmx(i, j, nz - ng) * 0.5d0

      fzmy_bnd(i, j, 1) = fzmy(i, j, ng) * 0.5d0
      fzmy_bnd(i, j, 2) = fzmy(i, j, nz - ng) * 0.5d0

      fzmz_bnd(i, j, 1) = fzmz(i, j, ng) * 0.5d0
      fzmz_bnd(i, j, 2) = fzmz(i, j, nz - ng) * 0.5d0

      fzbx_bnd(i, j, 1) = fzbx(i, j, ng) * 0.5d0
      fzbx_bnd(i, j, 2) = fzbx(i, j, nz - ng) * 0.5d0

      fzby_bnd(i, j, 1) = fzby(i, j, ng) * 0.5d0
      fzby_bnd(i, j, 2) = fzby(i, j, nz - ng) * 0.5d0

      fzbz_bnd(i, j, 1) = fzbz(i, j, ng) * 0.5d0
      fzbz_bnd(i, j, 2) = fzbz(i, j, nz - ng) * 0.5d0

      fzen_bnd(i, j, 1) = fzen(i, j, ng) * 0.5d0
      fzen_bnd(i, j, 2) = fzen(i, j, nz - ng) * 0.5d0

      fzpo_bnd(i, j, 1) = fzpo(i, j, ng) * 0.5d0
      fzpo_bnd(i, j, 2) = fzpo(i, j, nz - ng) * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  second step of 2nd order Runge Kutta method
!----------------------------------------------------------------------------------------------------------------------------------|
   call weno_lx(fxro, fxmx, fxmy, fxmz, fxbx, fxby, fxbz, fxen, fxpo,                                                              &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre, nx, ny, nz, ch2)
   call weno_ly(fyro, fymy, fymx, fymz, fyby, fybx, fybz, fyen, fypo,                                                              &
      ro_pre, my_pre, mx_pre, mz_pre, by_pre, bx_pre, bz_pre, en_pre, po_pre, nx, ny, nz, ch2)
   call weno_lz(fzro, fzmz, fzmy, fzmx, fzbz, fzby, fzbx, fzen, fzpo,                                                              &
      ro_pre, mz_pre, my_pre, mx_pre, bz_pre, by_pre, bx_pre, en_pre, po_pre, nx, ny, nz, ch2)

   do k = 2, nz
      km1 = k - 1
      do j = 2, ny
         jm1 = j - 1
         do i = 2, nx
            im1 = i - 1
            ro(i, j, k) = 0.5d0 * (ro(i, j, k) + ro_pre(i, j, k) - dt_dx * (fxro(i, j, k) - fxro(im1, j, k)) -                     &
                                                                   dt_dy * (fyro(i, j, k) - fyro(i, jm1, k)) -                     &
                                                                   dt_dz * (fzro(i, j, k) - fzro(i, j, km1)))
            mx(i, j, k) = 0.5d0 * (mx(i, j, k) + mx_pre(i, j, k) - dt_dx * (fxmx(i, j, k) - fxmx(im1, j, k)) -                     &
                                                                   dt_dy * (fymx(i, j, k) - fymx(i, jm1, k)) -                     &
                                                                   dt_dz * (fzmx(i, j, k) - fzmx(i, j, km1)))
            my(i, j, k) = 0.5d0 * (my(i, j, k) + my_pre(i, j, k) - dt_dx * (fxmy(i, j, k) - fxmy(im1, j, k)) -                     &
                                                                   dt_dy * (fymy(i, j, k) - fymy(i, jm1, k)) -                     &
                                                                   dt_dz * (fzmy(i, j, k) - fzmy(i, j, km1)))
            mz(i, j, k) = 0.5d0 * (mz(i, j, k) + mz_pre(i, j, k) - dt_dx * (fxmz(i, j, k) - fxmz(im1, j, k)) -                     &
                                                                   dt_dy * (fymz(i, j, k) - fymz(i, jm1, k)) -                     &
                                                                   dt_dz * (fzmz(i, j, k) - fzmz(i, j, km1)))
            bx(i, j, k) = 0.5d0 * (bx(i, j, k) + bx_pre(i, j, k) - dt_dx * (fxbx(i, j, k) - fxbx(im1, j, k)) -                     &
                                                                   dt_dy * (fybx(i, j, k) - fybx(i, jm1, k)) -                     &
                                                                   dt_dz * (fzbx(i, j, k) - fzbx(i, j, km1)))
            by(i, j, k) = 0.5d0 * (by(i, j, k) + by_pre(i, j, k) - dt_dx * (fxby(i, j, k) - fxby(im1, j, k)) -                     &
                                                                   dt_dy * (fyby(i, j, k) - fyby(i, jm1, k)) -                     &
                                                                   dt_dz * (fzby(i, j, k) - fzby(i, j, km1)))
            bz(i, j, k) = 0.5d0 * (bz(i, j, k) + bz_pre(i, j, k) - dt_dx * (fxbz(i, j, k) - fxbz(im1, j, k)) -                     &
                                                                   dt_dy * (fybz(i, j, k) - fybz(i, jm1, k)) -                     &
                                                                   dt_dz * (fzbz(i, j, k) - fzbz(i, j, km1)))
            en(i, j, k) = 0.5d0 * (en(i, j, k) + en_pre(i, j, k) - dt_dx * (fxen(i, j, k) - fxen(im1, j, k)) -                     &
                                                                   dt_dy * (fyen(i, j, k) - fyen(i, jm1, k)) -                     &
                                                                   dt_dz * (fzen(i, j, k) - fzen(i, j, km1)))
            po(i, j, k) = 0.5d0 * (po(i, j, k) + po_pre(i, j, k) - dt_dx * (fxpo(i, j, k) - fxpo(im1, j, k)) -                     &
                                                                   dt_dy * (fypo(i, j, k) - fypo(i, jm1, k)) -                     &
                                                                   dt_dz * (fzpo(i, j, k) - fzpo(i, j, km1)))
         enddo
      enddo
   enddo

   call resistivity(bx_resis, by_resis, bz_resis, en_resis,                                                                        &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, x, y, z, dx, dy, dz, nx, ny, nz)

   call source_term(mx, my, mz, bx, by, bz, en, po,                                                                                &
      ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, po_pre, gx, gy, gz, bx_resis, by_resis, bz_resis, en_resis,          &
      ch2, cp2, dx, dy, dz, dt / 2.0d0, nx, ny, nz)

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value in x direction
!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
      fxro_bnd(j, k, 1) = fxro_bnd(j, k, 1) + fxro(ng, j, k) * 0.5d0
      fxro_bnd(j, k, 2) = fxro_bnd(j, k, 2) + fxro(nx - ng, j, k) * 0.5d0

      fxmx_bnd(j, k, 1) = fxmx_bnd(j, k, 1) + fxmx(ng, j, k) * 0.5d0
      fxmx_bnd(j, k, 2) = fxmx_bnd(j, k, 2) + fxmx(nx - ng, j, k) * 0.5d0

      fxmy_bnd(j, k, 1) = fxmy_bnd(j, k, 1) + fxmy(ng, j, k) * 0.5d0
      fxmy_bnd(j, k, 2) = fxmy_bnd(j, k, 2) + fxmy(nx - ng, j, k) * 0.5d0

      fxmz_bnd(j, k, 1) = fxmz_bnd(j, k, 1) + fxmz(ng, j, k) * 0.5d0
      fxmz_bnd(j, k, 2) = fxmz_bnd(j, k, 2) + fxmz(nx - ng, j, k) * 0.5d0

      fxbx_bnd(j, k, 1) = fxbx_bnd(j, k, 1) + fxbx(ng, j, k) * 0.5d0
      fxbx_bnd(j, k, 2) = fxbx_bnd(j, k, 2) + fxbx(nx - ng, j, k) * 0.5d0

      fxby_bnd(j, k, 1) = fxby_bnd(j, k, 1) + fxby(ng, j, k) * 0.5d0
      fxby_bnd(j, k, 2) = fxby_bnd(j, k, 2) + fxby(nx - ng, j, k) * 0.5d0

      fxbz_bnd(j, k, 1) = fxbz_bnd(j, k, 1) + fxbz(ng, j, k) * 0.5d0
      fxbz_bnd(j, k, 2) = fxbz_bnd(j, k, 2) + fxbz(nx - ng, j, k) * 0.5d0

      fxen_bnd(j, k, 1) = fxen_bnd(j, k, 1) + fxen(ng, j, k) * 0.5d0
      fxen_bnd(j, k, 2) = fxen_bnd(j, k, 2) + fxen(nx - ng, j, k) * 0.5d0

      fxpo_bnd(j, k, 1) = fxpo_bnd(j, k, 1) + fxpo(ng, j, k) * 0.5d0
      fxpo_bnd(j, k, 2) = fxpo_bnd(j, k, 2) + fxpo(nx - ng, j, k) * 0.5d0
   enddo
   enddo

   do k = 1, nz
   do i = 1, nx
      fyro_bnd(i, k, 1) = fyro_bnd(i, k, 1) + fyro(i, ng, k) * 0.5d0
      fyro_bnd(i, k, 2) = fyro_bnd(i, k, 2) + fyro(i, ny - ng, k) * 0.5d0

      fymx_bnd(i, k, 1) = fymx_bnd(i, k, 1) + fymx(i, ng, k) * 0.5d0
      fymx_bnd(i, k, 2) = fymx_bnd(i, k, 2) + fymx(i, ny - ng, k) * 0.5d0

      fymy_bnd(i, k, 1) = fymy_bnd(i, k, 1) + fymy(i, ng, k) * 0.5d0
      fymy_bnd(i, k, 2) = fymy_bnd(i, k, 2) + fymy(i, ny - ng, k) * 0.5d0

      fymz_bnd(i, k, 1) = fymz_bnd(i, k, 1) + fymz(i, ng, k) * 0.5d0
      fymz_bnd(i, k, 2) = fymz_bnd(i, k, 2) + fymz(i, ny - ng, k) * 0.5d0

      fybx_bnd(i, k, 1) = fybx_bnd(i, k, 1) + fybx(i, ng, k) * 0.5d0
      fybx_bnd(i, k, 2) = fybx_bnd(i, k, 2) + fybx(i, ny - ng, k) * 0.5d0

      fyby_bnd(i, k, 1) = fyby_bnd(i, k, 1) + fyby(i, ng, k) * 0.5d0
      fyby_bnd(i, k, 2) = fyby_bnd(i, k, 2) + fyby(i, ny - ng, k) * 0.5d0

      fybz_bnd(i, k, 1) = fybz_bnd(i, k, 1) + fybz(i, ng, k) * 0.5d0
      fybz_bnd(i, k, 2) = fybz_bnd(i, k, 2) + fybz(i, ny - ng, k) * 0.5d0

      fyen_bnd(i, k, 1) = fyen_bnd(i, k, 1) + fyen(i, ng, k) * 0.5d0
      fyen_bnd(i, k, 2) = fyen_bnd(i, k, 2) + fyen(i, ny - ng, k) * 0.5d0

      fypo_bnd(i, k, 1) = fypo_bnd(i, k, 1) + fypo(i, ng, k) * 0.5d0
      fypo_bnd(i, k, 2) = fypo_bnd(i, k, 2) + fypo(i, ny - ng, k) * 0.5d0
   enddo
   enddo

   do j = 1, ny
   do i = 1, nx
      fzro_bnd(i, j, 1) = fzro_bnd(i, j, 1) + fzro(i, j, ng) * 0.5d0
      fzro_bnd(i, j, 2) = fzro_bnd(i, j, 2) + fzro(i, j, nz - ng) * 0.5d0

      fzmx_bnd(i, j, 1) = fzmx_bnd(i, j, 1) + fzmx(i, j, ng) * 0.5d0
      fzmx_bnd(i, j, 2) = fzmx_bnd(i, j, 2) + fzmx(i, j, nz - ng) * 0.5d0

      fzmy_bnd(i, j, 1) = fzmy_bnd(i, j, 1) + fzmy(i, j, ng) * 0.5d0
      fzmy_bnd(i, j, 2) = fzmy_bnd(i, j, 2) + fzmy(i, j, nz - ng) * 0.5d0

      fzmz_bnd(i, j, 1) = fzmz_bnd(i, j, 1) + fzmz(i, j, ng) * 0.5d0
      fzmz_bnd(i, j, 2) = fzmz_bnd(i, j, 2) + fzmz(i, j, nz - ng) * 0.5d0

      fzbx_bnd(i, j, 1) = fzbx_bnd(i, j, 1) + fzbx(i, j, ng) * 0.5d0
      fzbx_bnd(i, j, 2) = fzbx_bnd(i, j, 2) + fzbx(i, j, nz - ng) * 0.5d0

      fzby_bnd(i, j, 1) = fzby_bnd(i, j, 1) + fzby(i, j, ng) * 0.5d0
      fzby_bnd(i, j, 2) = fzby_bnd(i, j, 2) + fzby(i, j, nz - ng) * 0.5d0

      fzbz_bnd(i, j, 1) = fzbz_bnd(i, j, 1) + fzbz(i, j, ng) * 0.5d0
      fzbz_bnd(i, j, 2) = fzbz_bnd(i, j, 2) + fzbz(i, j, nz - ng) * 0.5d0

      fzen_bnd(i, j, 1) = fzen_bnd(i, j, 1) + fzen(i, j, ng) * 0.5d0
      fzen_bnd(i, j, 2) = fzen_bnd(i, j, 2) + fzen(i, j, nz - ng) * 0.5d0

      fzpo_bnd(i, j, 1) = fzpo_bnd(i, j, 1) + fzpo(i, j, ng) * 0.5d0
      fzpo_bnd(i, j, 2) = fzpo_bnd(i, j, 2) + fzpo(i, j, nz - ng) * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine weno
